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Abstract. In this paper hyperbolic partial differential equations with random coefficients are 
discussed. We consider the challenging problem of ffux functions with coefficients modeled by spa- 
tiotemporal random fields. Those fields are given by correlated Gaussian random fields in space and 
Ornstein-Uhlenbeck processes in time. The resulting system of equations consists of a stochastic dif¬ 
ferential equation for each random parameter coupled to the hyperbolic conservation law. We define 
an appropriate solution concept in this setting and analyze errors and convergence of discretiza¬ 
tion methods. A novel discretization framework, based on Monte Carlo Finite Volume methods, is 
presented for the robust computation of moments of solutions to those random hyperbolic partial 
differential equations. We showcase the approach on two examples which appear in applications: The 
magnetic induction equation and linear acoustics, both with a spatiotemporal random background 
velocity field. 
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1. Introduction. Hyperbolic partial differential equations with random data 
have been an active research field over the last decades. In ample situations measure¬ 
ments are not accurate enough to allow an exact description of a physical phenomena 
by a deterministic model. To account for this, uncertainty is introduced in the ap¬ 
propriate parameters and the distribution of the (now stochastic) solution is studied. 
In this paper we consider linear hyperbolic PDEs with time and space dependent flux 
functions. Important examples of such equations include linear elasticity, and linear 
shallow water equations, as well as the linear acoustics and the magnetic induction 
equation considered in this article. 

1.1. Linear hyperbolic conservation/balance laws. Many important phys¬ 
ical phenomena can be modeled by first order linear hyperbolic systems. We consider 
a system of the general form 


( 1 ) 


9tU(x, i)+div- (A(x,t) U(x, i)) =S(x, i), 

Tirv n'\ — TT„^'v'l 


V D 7“ 


with suitable boundary conditions. Here U(x, t) denotes the vector of conserved 
quantities at a point x in the domain D at time t, and AU = (A^^U, • • • , A^^U), 
for d G N, are the flux functions in the Xi,--- ,x^^-direction, respectively, where 
^xr . ^ linear maps. The matrix A is assumed to be diagonalizable with 

real eigenvalues. Well-posedness of linear non-autonomous systems of conservation 
laws is studied for instance in [14, 15, 18, 26, 47]. For a linear advection equation, 
even in the case of variable coefficients, the characteristics never cross [26]. Looking 
at the linear advection equation Ut + {xu)x = 0 in one (spatial) dimension, one can 
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prove that all characteristic converge to the origin. It becomes clear that solutions 
for t ^ oo may become unbounded and contain Dirac delta functions. 

In this article we treat the case where the flux functions A(x, t,V) depend on 
m G N coefficients V = that are modeled as correlated random fields 

in space, and stochastic processes in time, defined on a probability space P). 

We assume from now on that the random fields are (^-almost surely) differentiable 
with respeet to x, [2]. The stochastic processes considered provide the coefficients 
implicitly, namely as the solution to a stochastic differential equation (SDE). This 
means that the system we are considering is given by 
( 2 ) 

r5tU(x,i) + div ■ ^A(x,t, V(x,t)) U(x,i)j = S(x,i), 

\atV(x,t) = |i(V(x,i))+<7(V(x,i)) i^(x,t), 

where the vectors /[x, and cr are arbitrary functions, and the vector i/ is a random 
funetion in spaee and time, often referred to as the noise term”. 

In general, the Itb-SDE for the parameter V in Equation (2) has no explicit 
(strong or weak) solution. Many stochastic schemes fall into the class of stochastic 
Runge-Kutta (SRK) methods. Weak approximations focus on the expectation of 
functionals of solutions, whereas strong approximations are concerned with pathwise 
solutions. Eor an overview on the theory of SDEs and numerical methods we refer to 
[24, 23, 36, 29, 41] and references therein. 

In general, there are no explicit solution formulas for deterministic variable- 
coefficient linear hyperbolic systems of the form (1), let alone for the stochastic PDE 
(2). Numerical methods are therefore widely used to approximate the solutions of 
those types of equations. Einite Difference, Einite Volume and Discontinuous Galerkin 
methods are popular approaches to obtain efficient time and space discretizations for 
the probiem at hand, see [26, 18] and references therein. 

1.2. Uncertainty quantification. In many appiications the parameters of the 
fiux functions in Equation (1) are determined by measurements. Then is, in fact, only 
statistical information available. Among other phenomena, scarcity of measurements 
of material properties or background velocity fields lead to uncertainty in the param¬ 
eters of the flux function. Given the statistical description of the parameters, it is 
of interest to efficiently quantify the resulting uneertainty in the solution to Equa¬ 
tion (1). An appropriate mathematical notion, together with a proof of existence 
and uniqueness, of random solutions for systems of hyperbolic conservation laws has 
recently been developed in e.g. [4, 42] (see also references therein). 

Efficient numerical methods for uncertainty quantification in the setting of partial 
differential equations has been intensively studied in recent years. A non-exhaustive 
list of literature on uncertainty quantification for hyperbolic conservation laws in¬ 
cludes [1, 6, 27, 28, 44, 38, 46, 16, 21, 43, 38] and references therein. Among the 
most popular techniques are stochastic Galerkin methods based on generalized poly¬ 
nomial ehaos expansions (gPC). These methods reduce the stochastic model to a 
(high-dimensional) deterministic one. This comes to the price that they are highly 
intrusive, such that existing numerical schemes for conservation laws cannot be used. 
A second class of methods for uncertainty quantification are stochastic collocation 
methods, which are non-intrusive and easier to parallelize than gPC based methods. 
Solutions to hyperbolic conservation laws, however, do not have the necessary regu¬ 
larity with respect to the stochastic variables, which in general diminishes the use of 
both gPC and collocation methods. There are a number of other techniques, namely 
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stochastic Finite Volume methods (see [30]), adaptive analysis of variance, proper 
generalized decomposition, and Fokker-Planck-Kolmogorov type techniques. The lat¬ 
ter can handle low parametric regularity, but assume that the ’’effective” number of 
stochastic dimensions is low and require often impractical complex representations of 
the input random fields. 

Here, we use Monte Carlo (MC) methods to quantify the uncertainty, which com¬ 
prises a class of non-intrusive methods well suited for problems with low parametric 
(stochastic) regularity. Monte Carlo methods rely on repeated sampling of the prob¬ 
ability/parameter space. For each sample the underlying, (then) deterministic, PDF 
is solved and the sample solutions are combined to obtain statistical information in 
the form of moments of the distribution. Monte Carlo methods are very robust with 
respect to the regularity of the solutions. However, the convergence rate of a (plain) 
MC method is limited to ^ with respect to the number of samples, which means that 
typically a large number of realizations of approximations of solutions to the underly¬ 
ing problem (1) has to be computed. In the Monte Carlo approximation of moments 
of solutions to stochastic partial differential equations the discretization error consists 
of a statistical error and a spatiotemporal error of the numerical scheme, as can be 
seen in Equation (31). For non-autonomous linear systems of conservation laws it has 
been shown in [34] , that the number of Monte Carlo samples M can be chosen 

(3) M = 

in order to equilibrate the statistical and spatiotemporal error of an underlying FV 
scheme with convergence rate o > 0. The computational complexity due to the slow 
convergence rate, can be lowered by Multilevel Monte Carlo (MLMC) methods, which 
have been proposed in [42] and related papers by the same authors [33, 31, 32]. For 
a result on the convergence and computational complexity of the multilevel Monte 
Carlo approximation for general Hilbert-space-valued random variables see [3]. The 
idea behind MLMC methods is to use a hierarchy of different levels of resolution 
of the underlying deterministic numerical solver. The level dependent numbers of 
MC samples are then chosen, in such a way that the total computational work is 
minimal while the sum of all error terms is asymptotically optimized. In the case 
of non-autonomous linear systems of conservation laws, see [34] for details. Optimal 
computational complexity is not the main focus of this paper. We remark that it is 
not particularly challenging to employ an MLMC method, as the problems considered 
here fall into the class of problems for which a general MLMC method is derived in [3] . 

Random (scalar) linear transport equations are discussed for instance in [37, 10, 7, 
39] and references therein. In [11] the authors present expressions for the distribution 
of the solution of a linear advection equation with a time-dependent velocity, given in 
terms of the probability density function of the underlying integral of the stochastic 
process. Numerical schemes are then introduced in [9, 12]. A stochastic collocation 
method for the wave equation is introduced in [35]. Uncertainty quantification of 
acoustic wave propagation in random heterogeneous layered media is presented in [34] . 
In [21] the linear advection equation with spatiotemporal coefficients is the subject of 
research. The authors develop numerical methods using polynomial chaos expansions 
to solve the advection equation with a transport velocity given by a Gaussian or a 
log-normal distribution. 

1.3. Aims and contributions of the paper. We extend the existing frame¬ 
work for random solutions to linear hyperbolic systems presented in, e.g., [42, 33, 31, 
32] and the references therein. Random coefficients of the numerical fluxes are mod- 
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eled as Gaussian random fields in space and Ornstein-Uhlenbeck processes in time. 
Thus, the system of equations consists of a stochastic differential equation (SDE) for 
each random parameter coupled to the conservation law. For that purpose, 

• we provide the necessary solution concepts for weak random solutions for lin¬ 
ear conservation laws with time-dependent random field coefficients, together 
with a well-posedness result (existence, uniqueness and dependence on initial 
data). 

• we describe algorithms for fast generation of such spatiotemporal random 
fields. Discretizations of the SDEs with the appropriate (weak) order are 
presented. 

• we present two applications: the magnetic induction equation and linear 
acoustics, both with background velocity fields modeled by spatiotemporal 
random fields. 

The simulation results are based on a flexible, thread-parallel algorithm for the un¬ 
certainty quantification of linear conservation laws. 

The remainder of the paper is organized as follows. In Section 2 , we provide the 
necessary theoretical background for time-dependent random fields, as well as ran¬ 
dom (linear) hyperbolic conservation laws. We present a Monte Carlo Finite Volume 
framework for approximating moments of the random solutions in Section 3. This is 
followed by Section 4, showcasing efficiency and robustness on realistic test cases, and 
some conclusions in Section 5. 

2. Stochastic linear hyperbolic conservation/balance laws with spa¬ 
tiotemporal random parameters. This article treats the case of random flux 
functions with coefficients that are modeled as spatiotemporal random fields. In order 
to rigorously define uncertainty of the parameters of the flux function, we start by 
recapitulating the necessary background from probability theory. 

2.1. Spatiotemporal random fields. To this end, let be a measurable 

space with elementary events uj e ft endowed with a a-algebra A. Then, the measure 
space (f], M, P) is called a probability space^ if P : D ^ [ 0 , 1 ] is a a-additive set function 
such that P(f]) = 1. 

Definition 1 (Gaussian random field (GRF)). A random field g = { 5 f(x),x G 
D} (also written as g(K^uj)) is a set of real-valued random variables defined on a 
probability space (f2,M, P). The field g(K^uj) is called Gaussian, if the vector of ran¬ 
dom variables follows the multivariate Gaussian distribution for any x G D, i.e., 
g ^ where J\f is the normal distribution with mean {i and covariance ma¬ 

trix C(x, y). Any Gaussian random field is completely defined by its second-order 
statistics. 

We note that C is a nonnegative, semi-definite, symmetric function. Bochner’s 
theorem [5] states that C is the Fourier transform of a positive measure pic on 
If we assume that pic has a Lebesgue density 7 which is even and positive, we can 
construct a GRF by 

(4) 5(x) = 

where T denotes the d-dimensional Fourier transform with inverse ^ and W is 
a centered Gaussian family W = {IT(x),x G with covariance E[IT(x)IT(y)] = 
(5(x-y),x,y G (see [25]). 

Looking at the time domain, a standard Brownian motion or Wiener process 
B = {B{t)A ^ T < c>c, defined on a probability space (r^,M,P), is a 
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(a) Ornstein-Uhlenbeck process in time. (b) Correlated Gaussian random field in space. 


Fig. 1. Sample solutions of random processes in time and random fields in space. 


continuous stochastic process which starts in zero P-a.s. and has independent and 
normally distributed increments, i.e., Bt — Bg ^ — s). 

Definition 2 (Ornstein-Uhlenbeck (OU) process). Let B = {B{t),t e [0,T]) 
be a standard Brownian motion. For /i, 6>, a G M, > 0 and cr > 0, the Ornstein- 
Uhlenbeek proeess is given as the solution of the stoehastie differential equation 

da{t) = — a{t))dt F (jdB{t)^ 

a(0) = ao, 

In general the initial eondition may be random as well. 

The OU processes is mean-reverting: Starting at a value ao, over time, the process 
tends to drift towards its long-term mean ja. See Figure 1 for some sample solutions. 

Remark 3. For every t G [0,T] the random variable a{t) is normally distributed 
with mean and varianee 


E(a(t)) =/i + (ao -/i)e 

( 6 ) 

V(a(t)) = ^(l-e-^^^). 

This ean easily be shown by using Ito’s formula with the funetion f{t,x) = e^^x, and 
eonsidering the dynamies of f{t^a{t)). Then, the stoehastie differential equation (5) 
has the following solution 

(7) a(t) = /i + e“^^(ao — /i) + cr f dB{s). 

Jo 

From this form we ean direetly deduee the expeetation ofa{t), the varianee is derived 
by using the ltd isometry. 

Remark 4. The realizations of an Ornstein-Uhlenbeek proeess are eontinuous 
and nowhere differentiable with probability 1. 

Definition 5 (Spatiotemporal random field). For all t G M+ letG{t) = {G(x,t),x g| 
D} be a Gaussian random field with eovarianee G and mean 0. Further, a time- 
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dependent random field is defined as the solution Z of the following SDE (ef. Defini¬ 
tion 2) 


(iZ(x, t) = 6>(/i(x) — Z(x,t))6it + cr(iG(x, t), 

Z(x,0) = Zo(x), 

where ji is a eontinuously differentiable funetion in L^{D), and O^a are positive real 
parameters. 

The solution Z : x ^ M has the following properties. 

• For a fixed time i G M+; {Z(x, t),x G D} is a real-valued Gaussian random 
field. 

• For eaeh point x G (Z(x, t),t G M+) is an Ornstein-Uhlenbeek proeess, 
i.e., a mean-reverting proeess. 

2.2. Random conservation laws. Equipped with a probability space (1^, T., P) 
we can incorporate uncertainties in the flux of Equation (1) by considering the equa¬ 
tion 


(9) 


dtU(x, t,cc;) + div • (A(x, t,cc;)U(x, t,cj)) = S(x, 

U(x,0,u;) = Uo(x,u;), 


X G C > 0. 


We are interested in cases, where some or all of the coefficients of A are modeled 
as a spatiotemporal random field according to Definition 5. In order for this to 
make sense, we follow [42, 34], but extend the solution concept and definition of 
hyperbolicity to be time-dependent where necessary. 

Definition 6 (Hyperbolicity). Forw in the unit sphere let A’^(x, t,ci;) = 
uj) be the eonvex eombinations of the direetional random matriees 
. Consider the eigen-deeomposition 

= Q^(x, t,co’)A^(x, t,co’)Q^(x, 

where is a diagonal matrix eonsisting of the eigenvalues {Xfi ^ 1 <r <d) of , 
and eontains the eorresponding eigenveetors as eolumns. The random linear sys¬ 
tem of eonservation laws (9) is P-a.5. hyperbolie if all eigenvalues of A^ are real 
¥-a.s. for all (x, t) G x M+. In addition, for every finite time horizon T < oc there 
exists a K{uj) < oc sueh that 

(10) sup ||Q^(x,t,co’)||||Q^(x, t,cc;)“^|| <K{uj), P-a.5.. 

xG-D,tG[0,T],wG§'^“^ 

In addition, we require the expected wave speeds to be finite, i.e., 

(11) A := max sup |E [A)^(x, t, •)] | < oc, VtG[0 ,T]. 

l<r<d xGT),wG§^“^ 

We consider a measurable mapping 


(12) U : (D, A) (E, S(E)), cj ^ U(x, t, uj), 

where B{V) is a Borel cr-algebra of the function space V. Then, we define the concept 
of a weak solution as follows. 
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Definition 7 (Pathwise weak solution). A random field \J, with values in 
C([0, T], L^(D))^ is ealled a (pathwise) weak solution to the stoehastie eonservation 
law (9) on D X [0, T], with D = and a finite time horizon T < oo, if it is ¥-a.s. a 
weak solution. This means that it satisfies the variational formulation 

( 13 ) 

dxidt / Uo • = / S-(pdxidt, 

JR^X[0,T1 \ J JRd J^a 


/R^x [0,T] 


for all test funetions (p G x [0,T]) for F-a.e. uj G Q. 

Denote L^(D) = L^{D)^ and W’^’^(D) = {D)^. We have the following 

result regarding existence and uniqueness of a solution. 

Theorem 8 (Well-posedness of stochastic linear hyperbolic conservation laws.). 
Consider the linear eonservation law (2), and assume the following holds: 

• the stoehastie differential equation for the eoeffieients V (see Equation (20)) 
admits a unique solution (in the sense of [22, Chapter 5.2]) 

• the system is hyperbolie aeeording to Definition 6, with (pathwise) eonstant 
K = ||A(cj)||Lfc('Q^]^) < OC; 

• the moments of the initial data, the souree and the flux are hounded in the 
following sense: there exist non-negative r^^rs^VA G NU {0, oc} sueh that 

(14) Uo G L^(fl,W^°’^), S G G L^(fl,W^^’^), 

where Tt is the sample spaee of the probability spaee. 

• eaeh random field is stoehastieally independent o/Uq and S. 

Then, for T < oc, the system (2) admits a unique pathwise weak solution. Moreover, 
for all t G [0,T]^ we have the following estimates 


(15) 


||U(-5^5^)||l2(d) < K{iu) (||Uo(.,u;)||l2(^) + ^||S(.)||l2(t))) ^ ^-a.s., 

l|U||Lfc(Q,C([0,T],L2(D)) < A (||Uo||Lfc(Q,L2(D)) + ^||S||Lfc(0,L2(D))) 5 


Proof. The proof is an immediate consequence from [42, Theorem 1], if one con¬ 
siders the following modifications to address the time-dependence of the flux function: 

1. Using the time-dependent version of the classical existence and uniqueness 

results summarized in [42, Theorem 1], one can show that the random field 
given by the uj U(-, •, cj) is well defined and that is a weak solution 

P-a.s.. 

2. To show that the maps uj ^ \Jf,t,uj) are measurable for all t G [0,T] P-a.s. 
one uses the fact that LA{D) is a separable Hilbert space and the stability 
estimate of the classical theorem for the deterministic (pathwise) solution. 

3. The first estimate in (15) may again be derived from the time-dependent 
version of the classical case (summarized in [42, Theorem 1]), and the sec¬ 
ond follows from the first, using the assumption that random fields, initial 
conditions and sources are independent. 

Incorporating the above into the proof of [42, Theorem 1], the assertion follows im¬ 
mediately. □ 

If Uo,S G I/^(U,L^(D)) and K G Theorem 8 ensures the existence of 

moments of order k of the random weak solution. 

In general, there are no explicit solution formulas. For the special case of the 
scalar linear transport equation with a time-dependent coefficient (transport driven 


7 


by the Ornstein-Uhlenbeck process) however, we can derive the distribution of the 
solution in closed form. This distribution will then be used in an example presented 
in Section 4.1 to verify our Monte Carlo Finite Volume method. 

Theorem 9. Consider the sealar transport equation 


(16) 


u{x^t^uj)t + {a{t^uj)u{x^t^uj))x = 0 

0 ) = uq{x) 


with eoeffieient a{t^-) given by the Ornstein-Uhlenbeek proeess (5). The moments of 
the solution to Equation (16) exist and are given by 


(17) E{u{x,t)) = j Uo{x-y)fA(a^,ii){y)dy = {fA(a^,0)*Uo){x-tfl). 

Here, denotes eonvolution and the probability density funetion /a is given by 


fA{a‘^ r— 

y 


-e 26-2 


^TTCr^ 

3“^^ — — I) and transportation speed 

Remark 10. Higher moments of the solution may be ealeulated by 


with diffusion eoeffieient = p- ( 6 >t + 2 e — ^e 


(18) M^(R(x,t)) =E((^(x,t) -E(^(x,t)))"^). 

Proof The solution for a single realization (fixed u; G 11) of Equation (16) is given 
by uo{x — Jq a{s,uj) ds). We start by calculating the first moment of this expression, 

i.e. 


E{uo{x — / a{s)ds)). 

Jo 

This means, we have to calculate the distribution of the time integral over a, i.e. the 
distribution of the stochastic process 

A{t) = f a{t) dt. 

Jo 

The process A is again a Gaussian process, i.e. A{t) ^ and therefore 

completely characterized by its mean and variance. Using Fubini’s theorem we have 
that 


(19) 

C C — 1 

E{A{t)) = / E{a{s)) ds = / la P (oq — la) ds = ffi — {oq —/a) - - - =: fi. 

Jo Jo ^ 

We express the variance of A via the covariance of A with itself 

Y{A(t)) = CoY{A{t), A{t)) = E((A(i) - E{A{t))){A{t) - E(A(t)))). 
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Using A{t) —E{A{t)) = a Jq/qG dB{u)ds (combine Equations (7) and (19)) 


this yields 



Y{A{t))=E{a / e-^^^-^UB{u)dsa / e-^^^-^UB{v) dr) 


Jo Jo Jo Jo 



using Fubini’s theorem. For a Brownian motion 5, it is known that 



Therefore, we have 



This gives us the variance of A{t) depending on the variables 0 and a. 

Therefore, the expectation of the solution to Equation (16) is given by 



where /a is the normal density function with parameters fi and given by fA{y) = 



1 


□ 




Remark 11 . For the limit 0^0, we reeover the eorresponding result for a pure 


Brownian motion proeess (i.e. a{t) = (jB(t)), where jl = /i and = cr^y. This ean 
he shown by a Taylor expansion as 


Y{A{t)) 



A similar Taylor expansion shows the result for E{A{t)). 

3. Monte Carlo Finite Volume methods. For the approximation of the (mo¬ 
ments of the) solution to partial differential equations with random coefficients we 
have to discretize in space and time, as well as in the “stochastic domain”. We use a 
Monte Carlo method for the approximation of moments of the random solution. This 
means that we have to approximate the (deterministic) solution for each realization 
uj of Equation (9), i.e.. 


( 20 ) 



dA (x, t) = n{V (x, t)) + cr(V (x, t)) i/(x, t) 
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where each component of V is given by the spatiotemporal random field (8). As 
mentioned in the Introduction, a Multilevel MC method could be used to achieve 
optimal computational complexity. This is, however, not the goal of this article. 

Before we start with a brief recapitulation of Finite Volume methods, we introduce 
some useful notation. Let the computational domain be a bounded axiparallel domain, 
i.e., D = Dx X Dy X Dz. A uniform axiparallel mesh of the domain D consists of 
identical cells ^OTl<i<I,l<j<J,l<k<K and I^J^K < oc, with 

edge lengths Ax, A^, A^ in the x-, y- and z-direction, respectively. The cell centers 
are denoted by and values at the cell interfaces by ^i-i/ 2 : ^j- 1/25 ^/c-i/ 2 - For 

a function /(x, t) : !d ^ M we set /T ^ where for n = 1 ,..., V, is 

the n-th time step. 

3.1. Finite volume methods. A Finite Volume (FV) scheme is obtained by 
integrating Equation (1) over a cell, or control volume, and over a time in¬ 
terval + At'^. Denoting cell averages by = 

1 ^.^, ^ U(x, t) dx, a fully discrete flux-differencing method has the form 


At^ 


T jn+l ^ T jn _ / p7 

( 21 ) 

_ f^n 


_ rpn \ 

-|-l/2,jf,/c i—lj2,j,k) 


Af^ 

/^n \ _ ( TTn _ TT 

—1/2,/c/ \-^i,j,k-\-l/2 -^ijjfj/c —1/2 


- 1 / 2)5 


where F^G^H are the fluxes in x-, y- and z-direction respectively. The fluxes j k 

approximate the integral 


I r ryj+1/2 rZk+1/2 


The approximations for the fluxes in the other directions are equivalently defined. 
Using the flux-differencing form, FV methods are typically based on the reconstruct- 
evolve-aver age (REA) algorithm. This algorithm consists of the following steps per¬ 
formed for each time step: Eirst, one reconstructs the cell averages by a piecewise 
polynomial function. Then, one evolves the hyperbolic equation by defining the nu¬ 
merical fluxes (22). Einally, the new cell averages are computed. The numerical 
fluxes are based on solutions of local Riemann problems at each cell interface. High 
order accuracy in space is achieved by using TVD limiters, or (W)ENO schemes 
[19, 26, 45, 20, 40], and in time by strong stability-preserving (SSP) Runge-Kutta 
methods, see e.g., [17]. The CEL condition dictates a limit on the time step size At^ 
for the resulting explicit EV schemes. Let A^, A^, A^ be the eigenvalues of Equation (1) 
in X-, y-, and z-direction, then the time step size AV^ needs to satisfy 


(23) 


Ar^ Ar^ Ar^ 1 

Ax Ay Az 2 


where A^ = maxmax|AMx^ ^ t^)| is the largest absolute eigenvalue in x-direction, 

i,j,k p 

and similar in y-, and z-directions. 

3.2. Realizations of spatiotemporal random fields. In order to approxi¬ 
mate moments of the solution to Equation (9), we need to create realizations of the 
spatiotemporal random fields given in Definition 5. Therefore, we describe here one 
way to discretize the SDE (8) in time and space. 
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3.2.1. Discretization in time. There is a wide variety of numerical methods 
for SDEs. Methods are classified mainly according to their strong (ps) and weak {p^) 
orders of convergence. The order of weak convergence is typically higher than the order 
of strong convergence of a scheme. For example the Euler-Maruyama approximation 
has convergence orders (ps^Pw) = (1/2,1), and the Milstein method has convergence 
orders {ps^Pw) = (1? !)• Foi* a classification of higher order SRK schemes see, e.g., [ 8 ]. 

For the Ornstein-Uhlenbeck process, given in Equation ( 8 ), the Milstein method 
is equivalent to the Euler-Mar uyama method, since cr(V) = a is independent of the 
process. The Milstein method is given by 

(24) Z^+^(x) - Z\x.) = hO (/i(x) - Z\x.)) + aVhG\x.), 

where Z\^ = Z{'K^t^) for all x G D and G^(x) = is a Gaussian random field 

in X at the discrete times tK 

Depending on the parameters, the Ornstein-Uhlenbeck process, given in Equa¬ 
tion ( 8 ), is stiff and appropriate implicit schemes have to be used. For instance, the 
implicit Milstein method is given by 


(25) 


Z^+^(x) = ^^^(x) + h6>/i(x) + cr\/hG^(x)^ / (1 + hO ), 

=/i(x). 


In general, higher order SRK schemes become increasingly involved, for the 
Ornstein-Uhlenbeck process, however, those schemes simplify considerably, as the 
function cr(V) = cr is constant. 

3.2.2. Discretization in space. The discretization in Equation (25) is only 
semi-discrete as it is continuous in space. For a fully discrete approximation of a real¬ 
ization of Equation ( 8 ) we need to provide an algorithm that approximates realizations 
of the Gaussian random fields G^(x) on a Gartesian grid over D C at each time 
step. To this end, we use the approach described in [25], which provides the following 
algorithm based on the representation formula (4). For simplicity we only present the 
periodic case. Let be the discrete Fourier transform with inverse and let Zij^k 
be random samples from a normal distribution for all 0 < i < /, 0 < j < J^O < k < K, 
i.e., Zij^k ^ -^(0, |Ac|“^), where |Ac| is the volume of the cells. Then, a Gaussian 
random field is given by 

( 2 b) {^i,j,k} — ^ ^ , 

where the covariance is given by the Fourier transform of the symmetric, positive 
function 7 . A typical family of functions for the Lebesgue density 7 is given by 

(27) = (1 + l|Pi,i,fc|l 2 )“'> A:,; e N, > 1, 

where the points in the Fourier domain are given by Pi,j,/c = (i — //2, j — J/2, /c — 
Kl2ff. Another possibility would be to employ an exponential covariance function 
given by 7 i,j,/e = exp(—||p^^j^/c|| 2 /'r’), with correlation length v. 

3.3. Monte Carlo approach. We employ a Monte Garlo Finite Volume method 
(MG FV method) in order to approximate the moments of the stochastic PDE ( 2 ) 
where particular coefficients are modeled as spatiotemporal random fields given in 
Equation ( 8 ), such that the system is hyperbolic according to Definition 6 . As before, 
we consider a bounded axiparallel domain D together with a mesh T consisting of 
identical cells Gij^k- 
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The step size Af^ of the explicit PDE solver is given by the CFL condition (23) 
and depends on the eigenvalues which are influenced by the stochastic parameters V. 
On the other hand, the interval length h of the discretization of the SDE is constant 
for each realization. To achieve an optimal convergence rate, the approximation error 
of the SDE’s and the FV metlmd should be of the same order. If a bound of the 
maximum expected eigenvalue A, as defined in Equation (11), is available, one can 
choose 

At 

(28) h = c^, 

A 

where c > 0 is a constant. We require the discrete times of the FV scheme to be 
a subset of the discrete times of SDE approximation, i.e., {t^, n = 1,. .., A^} C {t^ = 
/•h, / = 1,...,L}. 

The MC FV algorithm consists of three main steps. 

1. Generate M realizations 1 <m < M) for each parameter that 

is modeled as a spatiotemporal random field (8): 

• Choose an appropriate interval h = T/L for the approximation of the 
SDE. 

• Generate (L — 1) -M (times the number of parameters) Gaussian random 
fields G\ j ^ on the given mesh T according to Equation (26). 

• Use an approximation scheme with the same (weak) order as the FV 
method. 

2. For each generated realization, the deterministic problem of the underlying 
hyperbolic conservation law is solved on the given mesh T. In the under¬ 
lying Riemann problem, the random field - is assumed to be piecewise 
constant on the time intervals [/h, (/ + 1)^]- For the sample we denote by 
u(x,r, ujrn) the exact pathwise solution at time T, and by U 7 -(x, T, u;^) the 
numerical approximation. 

3. The M approximations of the sample solutions, i.e., U 7 -(x, T, cj^) are used 
to approximate moments of the random solution field U(x, T, cj). 

One is particularly interested in the first two moments, i.e., the expectation E[U] 
and the variance V[U]. The sample mean of the approximate solutions is used to 
estimate the expectation given by 


1 

(29) Em[Vt]{^,T) = ^ F Ut(x,T,c^„). 

m=l 

Higher statistical moments of U, such as the variance Um[U 7 -](x, T), can be approx¬ 
imated similarly. The total approximation error of the expectation in the norm, 

i.e. 


(30) ^appr(T') — 


is bounded by the sum of the numerical approximation error 5num of the base method 
and the Monte Carlo error 5 mcm 

(31) ^appr(t) ^ ^num (t) + ^MCM(t), 
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where 


I 1 ^ \ 

T, — U(x^^j|-^/e, T, Co’ttt,)) , 

i,j,k m=l 

I 1 ^ 

l^*li.fe| ^ y I ^ ^ U(Xj T, Wto) — ]E[U] (Xjj-^J;, T) . 

i,j,k m=l 

Using the triangle inequality, it is trivial to show the relationship (31). If one uses 
the norm then equality holds. 

The estimate (31) shows that the approximation error is bounded by the domi¬ 
nating part of the sum of the numerical error and the Monte Carlo error. The Monte 
Carlo method converges with the rate 1/2 in the number of samples in mean square 
and is independent of the resolution of the grid, i.e. the size of Ax. On the other hand, 
a numerical base method of order o converges with 0{/S.x^) for each single realization, 
independent of the number of Monte Carlo samples. Therefore, equation (31) suggests 
that our Monte Carlo method is most efficient if Snum — ^mcm- For the according 
sample numbers in a MLMC approach we refer to [3] . In [34] for instance it is shown 
that this can be achieved if the number of Monte Carlo samples is M = ©(Ax”^*^), 
where o is the order of the FV method. 

Lemma 12. Assume that, 

• the eonditions of Theorem 8 are fulfilled for k >2, 

• the underlying numerieal FV seheme eonverges (under grid refinement) to the 
weak solution of Equation (1) with rate o> 0, 

• the numerieal seheme for the SDE eonverges at the same rate o > 0. 

Then, the MC FV method estimates Em[V^t] deserihed in this seetion eonverge to 
the first moment of the solution E[U] in mean square sense, as Ax ^ 0, with M = 
0(Ax-2^). 

Froof The assertion follows directly from the triangle inequality and the con¬ 
vergence of the corresponding FV scheme, together with the specific choice for the 
number of samples 

E||E[U] -EM[Ur]||' < ||E[U] -E[Ur]||' + E|XM[Ur] -E[Ur]||' 

< E||U - Urf + ^Var(Ur) 

= o(aA°). 

Here we used the standard convergence properties of the sequence of Monte Carlo 
estimators G N). The norm || • || is the canonical norm for the (pathwise) 

solution U. □ 

There is no dependence between different samples and therefore the described 
algorithm is trivial to parallelize. Our implementation distributes the workload on as 
many CPU threads as are available in the computing environment, achieving (triv¬ 
ially) optimal parallelization-efficiency. 

4. Examples. We evaluate the proposed approach for uncertainty quantifica¬ 
tion of linear hyperbolic conservation laws on a suite of test cases. We start with 
a “degenerate” case, where an autonomous scalar linear transport is driven by a 
(space-mdependent) Ornstein-Uhlenbeck process. We present error and convergence 


(32) 


^num(^) — 


^mcm(^) = 
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Fig. 2. Expectation and standard deviation of the solution of the OU-driven linear transport (33). 


analysis, based on the existence of the explicit solution formula (17). Then, we present 
two realistic cases of linear systems of conservation laws in two dimensions, namely 
the equations for linear acoustics and the motion of magnetic fields (induction equa¬ 
tion). In both cases, we model the (background) velocity field in x-, and y-direction 
as a spatiotemporal random field, and we show hyperbolicity of the system. The 
simulations show the robustness of the approach and reveal interesting features of the 
moments of the solutions. 

4.1. Ornstein—Uhlenbeck process driven scalar linear advection. We 

start by considering the scalar linear stochastic conservation law with a parameter 
given by the Ornstein-Uhlenbeck process (5), i.e., 

Ut + {a{t)u)x = 0, 

u{x,0,uj) =uo{x), xeD = [0,l] 
da{t) = 9{ii — a{t))dt + adB{t), 
a(0) = ao, 

with periodic boundary conditions for u. The eigenvalue of the PDE is the matrix 
itself, i.e., and the normalized eigenvector is 1. The system is 

hyperbolic according to Definition 6, since ||(3|||| = 1 and, using Equation (6), 
we have that the expected maximum eigenvalue 

(34) |E[A(x, t, cj)] I = |E[a(t, cj)] | = |/i + (ao — | < oo, for all t G M+ 

is finite. 

We use the derived explicit solution formula from Theorem 9 for the moments for 
the solution to this equation to show convergence of the proposed MC EV method. 
To this end, we compute the first two moments of Equation (33) at time t = 1, 
started with an initial condition consisting of a discontinuity, i.e. uq{x) = l[i_i i^i], 
where 1a is the characteristic function of the subset A. Eurthermore, we choose tfie 
deterministic initial condition of the OU process to be a( 0 ) = — | and (/i, ^,cr) = 
(|,20, \ ). A few sample solutions for these parameters are plotted in Eigure 1 (a). 
We can see in Eigure 2 that the expectation E{u) at time t consists of the initial 
function uq transported with speed ft and smoothed out wave fronts in accordance 
with Theorem 9. The largest values of the variance are located around the (smoothed 
out) discontinuities. 
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Fig. 3. Relative errors (in %) at time t = 1 for the OU-proeess-driven linear transport Equa¬ 
tion (33). Both the first and seeond order MC FV method eonverge to the exact solution. 


The first order scheme MC FV scheme uses a standard upwind discretization of 
the deterministic problem, and the Milstein scheme for the OU-process. The second 
order scheme consists of a minmod fiux-limiter in space and a second order strong 
stability preserving (SSP) Runge-Kutta time-stepping for the deterministic problem, 
together with a (weak) second order stochastic Runge-Kutta scheme for the OU- 
process. For the discrete time interval of the OU-process we use ^ ^ — 2Ax. 

The number of Monte Carlo samples is chosen as M = For the described 

setup. Figure 3 shows the relative approximation error of the first two moments of 
the solution. Both the first and second order MC FV method converge to the exact 
solution. The convergence order for the first moment is 5 ~ 1, and for the second 
moment is s ~ 0.7. The second order MC FV method has a smaller error constant 
compared to the first order method. The full convergence order s = 2 for the second 
order scheme is not achieved, since the deterministic solution consists of discontinuous 
piecewise linear data. 

4.2. Linear Acoustics in 2 dimensions. Sound waves can be described using 
the vector of conserved variables U = (p, rt, n)^, where p is the pressure, u is the 
velocity in x-direction, and v in y-direction. Given a background density po G M+ 
and a bulk modulus of compressibility Kq G M+, the dynamics are governed by the 
following system of equations 


U 


:+(a(o)u) +(a(i)u) =0, 


(35) 


(x, t, Lo) WyKo 

A"" = A’"(x, i,a;) = I Wx/po UQ{x,t,u)) 0 


^v/Po 

Uq = UqWx + VQWy, 


Uo{x,t,uj). 


where mo(x, t, •) is the stochastic background velocity field in x-direction, and uo(x, t, •) 
in y-direction. Both uq and vq are spatiotemporal random fields as specified in Defi- 
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(a) Pressure. (b) Velocity in x-direction. (c) Velocity in y-direction. 


Fig. 4. A sample solution of the linear aeousties equation with a stoehastie haekground veloeity 
field. The solution shows that the haekground veloeity field distorts the propagating waves. 


nition 5, i.e., they are given as the solution of the following SDEs 

duo{-K,t) =0u{fiu{^) - uo{-K,t))dt ^ audGu{y^,t), 
dvo{^,t) = - t;o(x,+ cr^(iG^(x,t). 


Defining the sound speed as cq = y^Kof po^ the eigensystem of (35) is given by 


(37) 


A^3(x, t, Uj) = (x, t, Uj) T Co, (x, t, Uj) = (x, t, u;), 




( — PoCo 0 PoCo\ 

-Wy Wx 
Wy Wx Wy ) 


where the rows of consist of the right eigenvectors of The eigenvectors are 
deterministic and are in fact the same as in the deterministic case. Therefore, the 
bound in Equation (10) naturally holds. Eurthermore, we have the following estimates 


|E[V3(^’*’‘^)]I ^ 2co + |E[uo(x,i,w)]| + |E[vo(x,i,w)]|, 
|E[Ay(x,i,w)]| < |E['Uo(x,i,w)]| + |E[vo(x,i,w)]|. 


Eor U{)^vo defined as in (8) these expectations are given by 

E[uo(x, t, w)] = /i„(x) + (ao(x) - /i„(x))e“®*, 

E[t)o(x,i,a;)] = /i„(x) + (6o(x) - /i„(x))e“®*. 

Using this, it is easy to show the bound (11), and therefore that the system (35) is 
hyperbolic according to Definition 6. 

We test our approach for an initial condition given by 

(40) Uo(x) = (p(x), u(x), v(x))^ = (0,0,0)^. 


We point out that pressure values do not have to be positive, because the system (35) 
is derived by linearizing around a background state (po,'^o,co)- This means that 
p, 14, and V describe the perturbation relative to the background state. The domain is 
given by D = [0,1]^ with Neumann boundary conditions on the right, top and bottom 
boundary. On the left boundary we have an acoustic source given by 


(41) 


U(0,p,t) 


(0, sin(47rt), 0)^, if|p—||<0.05 
(0,0,0)^, otherwise. 
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(a) Mean of pressure. 
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(b) Variance of pressure. 



(c) Mean of velocity in y-direction. 


0.2 0.4 0.6 0.8 1.0 

(d) Variance of velocity in y-direction. 


Fig. 5. Statistic mean and variance of a propagating sound wave computed with the proposed 
approach at time t = 1.5. Waves are introduced on the left boundary and propagate radially through 
the domain, showing a symmetric structure, eonsisting of smooth eircular wave fronts for pressure. 
The varianees are observed where the mean of the solution changes sign. 


The covariance for the stochastic background velocities ixq, vq is given by the Lebesgue 
density (27) with q = 2 and I = 4. For the initial condition of Equation (8) we 
use Z(x, 0) = 0, for the mean /i(x) = 0, and the remaining parameters are set to 
0 = a = 1. 

As a numerical scheme we use an approximate Riemann solver of the HLL-type 
(see [26]), resolving the outermost waves. Let and Fp^ji denote the left and 

right state and flux respectively. Then the numerical flux is given by 


(42) 


F«LL(UL,Ufl) = 


Fl, if f ^ 

F*, if Sl < f < Sr, 
Fr, ifSR<^, 


where sr = sr = and F>^ is determined from conservation leading to 


( 43 ) 


SrFr — SrFr + SlSr{IJr — Ul) 
Sr — sl 


For first and second order MC FV methods, approximations of the random back- 
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First order MC FV scheme First order MC FV scheme 


Ax = Ay 

1/16 1/32 1/64 1/128 Ax = Ay 

1/16 1/32 1/64 1/128 

Em[p\ 
Em [u] 
Em ['^i 

99.7 % 97.2 %87.7 %64.3 % Vm[p] 

99.3 % 96.6 %87.3 %64.4 % Vm[u\ 

101.9 %99.3 %87.4 %62.8 % Vm [v] 

99.8 % 98.8 %96.0 %85.2 % 

99.9 % 99.4 %97.8 %91.0 % 
100.0 %99.8 %98.6 %91.9 % 

Second order MC EV scheme Second order MC EV scheme 

Em[p\ 
Em [u] 
Em ['^i 

98.8 % 90.9 %60.8 %26.7 % Vm[p] 

98.6 % 90.0 %61.3 %27.6 % Vm[u\ 

101.1 %91.6 %58.7 %24.9 % Vm [v] 

Table 1 

92.8 % 93.8 %72.5 %40.1 % 
97.2 % 93.2 %79.7 %61.0 % 
99.6 % 97.5 %86.6 %63.1 % 


(Self) convergence to the solution with the second order scheme with Ax = Ay = 1/256. 


ground velocity field (36) are based on the Milstein scheme, and a (weak) second order 
stochastic Runge-Kutta scheme, respectively. For the given parameters, the discrete 
time interval of the OU-process is h = ^ min{Ax, Ay}. 

Results of the deterministic FV simulation are given in Figure 4. For the given 
scenario, a sample solution is plotted at time t = 1.5 with a second order minmod 
scheme on a 256 x 256 mesh. We can see that the waves enter the domain at the left 
boundary and propagate radially through the domain. As expected, the waves are 
distorted due to the presence of the random background velocity field. 

For the stochastic MC FV simulation we compute the mean and variance of the 
solution with a scheme of order o using M = O ((1/Ax)“^^) Monte Carlo samples (see 
Equation (3)). The structure of the mean of the propagating waves shown in Figure 5 
resembles the structure of the waves seen in the deterministic simulation of one sample 
shown in Figure 4. The sound waves introduced on the left boundary travel radially 
through the domain, showing a symmetric structure, consisting of smooth circular 
wave fronts for pressure. The largest values of the variance of a conserved quantity 
are observed around sign changes of the mean of the solution of that quantity. Table 1 
shows self-convergence of the proposed MC FV method to a reference solution on a 
256 X 256 grid. 

4.3. Magnetic Induction Equation in 2 dimensions. The magnetic induc¬ 
tion equation describes the evolution of a magnetic field U for a given velocity field 
V. We use the symmetric form of the equations, see [13], given by 


(44) 


dtV + div(V 0 U - U 0 V) = -Vdiv(U), 
U(x, 0) = Uo(x), with div(Uo) = 0. 


The equations have the intrinsic constraint that the divergence of the magnetic field 
U is preserved in time, i.e., dt(divU) = 0. Therefore, the system (44) is analytically 
equivalent to the conservative form without any source term. We consider the equation 
in two dimensions, where the components of the velocity field V = (rt, v) are given by 
spatiotemporal random fields as defined in Equation (8), i.e., as the solution of the 
following SDEs 


(45) 


(iM(x, t) = 9u{fiu{^) - ^(x, t))dt + audGu{x, t), 
dv{x,t) = dy{iJ,y{x) — v{x,t))dt + aydGv{x.,t). 
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(a) Magnetic field in x-direction. 



4 

0 




(b) Magnetic field in y-direction. 



(c) Velocity field in x-direction. 



(d) Velocity field in y-direction. 


Fig. 6. Sample solution at time t = 0.75 for the magnetie induetion equation. The initial 
magnetie field has been adveeted and distorted as a result of the stoehastie veloeity field. 


The eigensystem of the symmetric system (44) is given by 

(46) A^ 2 (x,i,w) = V(x,i,w)-w, Q'^ = . 

It is easy to show that the system is hyperbolic according to Definition 6. The 
bound (10) trivially holds, as Q is the matrix identity. Furthermore, we have that 

(47) |E[A^ 2 (x,t,cj)]| < |E[^o(x,t,cj)]| + |E[t;o(x,t,cj)]|. 

Using (39), we can easily show the bound (11) on the eigenvalues. 

Numerical schemes approximating the solutions of the induction equation (44) 
have to address the divergence constraint. Here, we will use the ’’stable upwind 
scheme” presented in [13]. To test our approach, we consider the equation on the 
domain D = [—^, ^]^, with periodic boundary conditions, and a divergence free initial 
magnetic field, given by a potential function H, i.e., 

(48) = {dyA{x,y),dxA{x,y)), withH(x,^) = sin(27rx) sin(27r^) + ^ - x. 
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The mean of the background velocity V(x, see Equation (8), is defined as 

/ cos(27rx) + 2 sin(27r^) sin(27rx) + 2 cos(27r^) 

(49) /i(x, y) = I 1 H ^ , 1 H ^ 

We set the initial condition of the velocity field, i.e., of Equation (8), to equal the 
mean, i.e., Z{x,y) = fi{x,y), and the remaining parameters to 6> = 1, and a = 10. 
The covariance is given by the Lebesgue density (27) with q = 2 and 1 = 4. The 
approximations of the random background velocity field (45) are based on the Milstein 
scheme. Eor the given parameters, the discrete time interval of the OU-process is 
h = I min{Ax, Ay}. 

Results of the deterministic EV simulation are shown in Eigure 6, together with 
the sample vector field. The first order stable upwind scheme is used to approximate 
the solution at time t = 0.75 on a 256 x 256 mesh. We see that the initial magnetic 
field has been advected and distorted due to the space- and time-dependent velocity 
field. 

Eor the stochastic MC EV simulation we compute the mean and variance of the 
solution using M = 100 (^) Monte Carlo samples (see Equation (3)). The struc¬ 
ture of the mean of the propagating waves shown in Eigure 7 resembles the structure 
of the waves seen in the deterministic simulation of one sample shown in Eigure 6. 
The values of the variance exhibit an interesting structure, which will be analyzed 
in a forthcoming paper. The table in Eigure 7 shows the expected convergence rate 
for the first order scheme, for both the first and the second statistical moment. The 
second statistical moment has a larger error constant compared to the first moment. 

5. Conclusions and Outlook. Linear systems of hyperbolic conservation laws 
with random coefficients are considered. Those coefficients are modeled as time- 
dependent random fields, leading to a coupled system consisting of the conservation 
law and stochastic differential equations for each of the parameters. An appropriate 
solution concept is developed and a Monte Carlo based algorithm is presented to 
approximate statistical moments of the solution. Important examples are presented, 
namely linear acoustics and magnetic induction with random velocity field coefficients. 
The results reveal interesting structures in the moments of the solution. Error and 
convergence analysis validate the proposed method. 

In the future, we plan to establish a rigorous error analysis and formal convergence 
proofs. Eurthermore, we will increase efficiency of the approach by developing highly 
parallel Multi Grid MC EV methods, utilizing the power of CPUs on coarse grid levels 
and graphics processing units (GPUs) on fine grid levels. This will further facilitate 
simulations of more complex problems in three dimension, and problems where the 
time-dependent random fields are given by more complicated SDEs, or even given by 
stochastic partial differential equations (SPDEs). 
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(a) Mean of magnetic field in x-direction. 



(b) Variance of magnetic field in x-direction. 



(c) Mean of magnetic field in y-direction. 



(d) Variance of magnetic field in y-direction. 


Self convergence to solution with Ax = Ay = 1/256 


Ax = Ay 

1/16 

1/32 

1/64 

1/128 

Em[Uo] 

73.0 % 

60.0 % 

41.9 % 

19.6 % 

Em[Ui] 

79.9 % 

61.3 % 

41.1 % 

18.8 % 

<1 

II 

H 

<1 

1/16 

1/32 

1/64 

1/128 

Vm[Uo] 

98.8 % 

97.7 % 

89.9 % 

62.5 % 

Vm[Ui] 

98.9 % 

97.1 % 

88.8 % 

82.0 % 


Fig. 7. Results of the computation of mean and variance using the first order stable upwind 
scheme. The magnetic field shows interesting features. The first moments of the solution converge 
to the reference solution at the expected rate. 
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